      Subroutine new_phi(r,alpha,theta,dtheta,C,phiu,rho0,
     + 			 LogInit,LogRapid)
      !*****************************************************************
      include 'accur.for'
      include 'statd.for'
      Parameter(KMax=(IMax+1)*(JMax+1))
      Parameter(nnzero=10*KMax)
      Parameter(LICN=7*nnzero,LIRN=7*nnzero/2)
      Parameter(NIN=5,NOUT=1)
      Dimension r(0:IMax),alpha(0:IMax),theta(0:JMax),
     +          a_1(0:IMax),b_1(0:IMax),c_1(0:IMax),
     +          a_2(0:JMax),b_2(0:JMax),c_2(0:JMax),
     +          b(KMax),bim(KMax),bip(KMax),bjm(KMax),bjp(KMax),
     +          bK0(KMax),bKA(KMax),bKB(KMax),
     +          b0(0:JMax),
     +          bphim(KMax),bphip(KMax),
     +          ff(KMax)
      Dimension kbphIm(KMax),kbphIp(KMax)
      Dimension ir0(nnzero),ic0(nnzero)
      Dimension A(LICN),
     +          W(KMax),
     +          ICN(LICN),
     +          IDISP(10),
     +          IKEEP(KMax,5),
     +          IRN(LIRN),
     +          IW(KMax,8),
     +          IVECT(nnzero),
     +          JVECT(nnzero)
      Logical GROW,LBLOCK,ABORT(4)
      Logical LogInit,LogRapid
      Common /phi/phi(0:IMax,0:JMax)
      Common /delta/delta
      Common /pi/pi
      Common /IB/IB/JB/JB/I0/I0

      dr=one/IA
      dalpha=one/r(IA-1)-one/r(IA)
      KA=kk(IA,JMax)
      KB=kk(IB,JB)
      K0=kk(I0,JMax)
      do i=1,KMax
          b(i)=zero
          bim(i)=zero
          bip(i)=zero
          bjm(i)=zero
          bjp(i)=zero
          ff(i)=zero
      enddo
      do i=1,IA-1
          a_1(i)=r(i)**(one-three*delta)/(delta*dr)**2
     +           *(half*(r(i-1)+r(i)))**(one+delta)
          c_1(i)=r(i)**(one-three*delta)/(delta*dr)**2
     +           *(half*(r(i)+r(i+1)))**(one+delta)
          b_1(i)=-a_1(i)-c_1(i)
      enddo
      do l=1,IMax-IA
          i=IMax-l
          a_1(i)=alpha(l)**(one+three*delta)/(delta*dalpha)**2
     +           *(half*(alpha(l)+alpha(l+1)))**(one-delta)
          c_1(i)=alpha(l)**(one+three*delta)/(delta*dalpha)**2
     +           *(half*(alpha(l-1)+alpha(l)))**(one-delta)
          b_1(i)=-a_1(i)-c_1(i)
      enddo
      do j=0,JMax
          if((j.eq.0).or.(j.eq.JMax)) then
              if(j.eq.0) then
                  c_2(j)=two**2/dtheta**2
                  b_2(j)=-c_2(j)
              else
                  a_2(j)=two/dtheta**2
                  b_2(j)=-a_2(j)
              endif
          else
              a_2(j)=one/(sin(theta(j))*dtheta**2)
     +               *sin(half*(theta(j-1)+theta(j)))
              c_2(j)=one/(sin(theta(j))*dtheta**2)
     +               *sin(half*(theta(j)+theta(j+1)))
              b_2(j)=-a_2(j)-c_2(j)
          endif
      enddo
c     === aproximation of delta-operator ==>
c     === intermal space ==>
      do i=1,IMax-1
          do j=0,JMax
	      k  =kk(i,j) 
	      b(k)=b_1(i)+b_2(j)/r(i)**(2*delta)
	      bim(k)=a_1(i)
              bip(k)=c_1(i)
	      if(j.ne.0) then
                  bjm(k)=a_2(j)/r(i)**(2*delta)
              endif
	      if(j.ne.JMax) then 
		  bjp(k)=c_2(j)/r(i)**(2*delta)
              endif
          enddo
      enddo
c     <== intermal space ===
c     === r=0 ==>
      b(kk(0,JMax))=-6d0/dr**(2*delta)
      do j=0,JMax
          theta1=half
     +           *(theta(max(0,j-1))+theta(j))
          theta2=half
     +           *(theta(j)+theta(min(JMax,j+1)))
	  b0(j)=6d0/dr**(2*delta)*(cos(theta1)-cos(theta2))
      enddo
      do j=0,JMax-1
	  k=kk(0,j)
	  b(k)=one
	  bjp(k)=-one
      enddo
c     <== r=0 ===
c     === alpha=0 ==>
      do j=0,JMax
	  b(kk(IMax,j))=one
      enddo
c     <== alpha=0 ===
c     <== aproximation of delta-operator ===
      do i=0,IA 
          do j=0,JMax
	    k=kk(i,j)
	    if(k.ge.kk(0,JMax)) then
	      psi=Fpsi(degree(r(i),delta)*sin(theta(j)))
              rho=EOS(phiu*(C-phi(i,j)-psi),1)/rho0
              drho_dH=EOS(phiu*(C-phi(i,j)-psi),2)
              drho_dphiu=drho_dH/rho0*(C-phi(i,j)-psi)
              drho_dC=drho_dH/rho0*phiu
              drho_dphi=-drho_dH/rho0*phiu
              drho_dpsi=-drho_dH/rho0*phiu
              dphiu_dC=-EOS(rho0,3)/(C-phi(I0,JMax))**2
              dphiu_dphi0=EOS(rho0,3)/(C-phi(I0,JMax))**2
              dC_dphi0=(-EOS(zero,3)/EOS(rho0,3))
     +                 /(one-EOS(zero,3)/EOS(rho0,3))
              dC_dphiA=one
     +                 /(one-EOS(zero,3)/EOS(rho0,3))
              dC_dpsiA=one
     +                 /(one-EOS(zero,3)/EOS(rho0,3))
              dpsi_dphiA=-Fpsi0(degree(r(i),delta)*sin(theta(j)))
              dpsi_dphiB= Fpsi0(degree(r(i),delta)*sin(theta(j)))
              dpsi_dpsi0=phi(IB,JB)-phi(IA,JMax)
              dpsiA_dphiA=-Fpsi0(one)
              dpsiA_dphiB= Fpsi0(one)
              dpsiA_dpsi0=phi(IB,JB)-phi(IA,JMax)
              dpsiA_dpsi0A=dpsiA_dpsi0
              if(LogRapid) then
                  if(JB.eq.0) then
                      iomb=0
                  else
                      iomb=IB
                  endif
                  om_=degree(r(i),delta)*sin(theta(j))
                  im=degree(om_,one/delta)/dr
                  ip=im+1
                  omm=im*dr
                  omp=ip*dr
                  xm=(omp-degree(om_,one/delta))/dr
                  xp=(degree(om_,one/delta)-omm)/dr
                  phiom_=xm*phi(im,JMax)+xp*phi(ip,JMax)
                  dpsi0_dphi0=-one/(phi(IA,JMax)-phi(iomb,JMax))
                  dpsi0A_dphi0=-one/(phi(IA,JMax)-phi(iomb,JMax))
                  dpsi0_dphiA=-(phiom_-phi(I0,JMax))
     +                         /(phi(IA,JMax)-phi(iomb,JMax))**2
                  dpsi0A_dphiA=-(phi(IA,JMax)-phi(I0,JMax))
     +                          /(phi(IA,JMax)-phi(iomb,JMax))**2
     +                          +one/(phi(IA,JMax)-phi(iomb,JMax))
                  if(JB.eq.0) then
	              dpsi0_dphiB=zero
                      dpsi0_dphi0= dpsi0_dphi0
     +                            +(phiom_-phi(I0,JMax))
     +                             /(phi(IA,JMax)-phi(iomb,JMax))**2
	              dpsi0A_dphiB=zero
                      dpsi0A_dphi0= dpsi0A_dphi0
     +                             +(phi(IA,JMax)-phi(I0,JMax))
     +                              /(phi(IA,JMax)-phi(iomb,JMax))**2
                  else
	              dpsi0_dphiB=(phiom_-phi(I0,JMax))
     +                            /(phi(IA,JMax)-phi(iomb,JMax))**2
	              dpsi0A_dphiB=(phi(IA,JMax)-phi(I0,JMax))
     +                             /(phi(IA,JMax)-phi(iomb,JMax))**2
                  endif
	          dpsi0_dphim=xm/(phi(IA,JMax)-phi(iomb,JMax))
	          dpsi0_dphip=xp/(phi(IA,JMax)-phi(iomb,JMax))
              else
                  dpsi0_dphi0=zero
                  dpsi0A_dphi0=zero
                  dpsi0_dphiA=zero
                  dpsi0A_dphiA=zero
                  dpsi0_dphiB=zero
                  dpsi0A_dphiB=zero
	          dpsi0_dphim=zero
	          dpsi0_dphip=zero
              endif
              ff(k)=4d0*pi
     +              *( rho
     +                -drho_dphiu
     +                 *( dphiu_dphi0*phi(I0,JMax)
     +                   +dphiu_dC
     +                    *( dC_dphi0*phi(I0,JMax)
     +                      +dC_dphiA*phi(IA,JMax)
     +                      +dC_dpsiA
     +                       *( dpsiA_dphiA*phi(IA,JMax)
     +                         +dpsiA_dphiB*phi(IB,JB)
     +                         +dpsiA_dpsi0A
     +                          *( dpsi0A_dphi0*phi(I0,JMax)
     +                            +dpsi0A_dphiA*phi(IA,JMax)
     +                            +dpsi0A_dphiB*phi(IB,JB)
     +                           )
     +                        )
     +                     )
     +                  )
     +                -drho_dC
     +                 *( dC_dphi0*phi(I0,JMax)
     +                   +dC_dphiA*phi(IA,JMax)
     +                   +dC_dpsiA
     +                    *( dpsiA_dphiA*phi(IA,JMax)
     +                      +dpsiA_dphiB*phi(IB,JB)
     +                      +dpsiA_dpsi0A
     +                       *( dpsi0A_dphi0*phi(I0,JMax)
     +                         +dpsi0A_dphiA*phi(IA,JMax)
     +                         +dpsi0A_dphiB*phi(IB,JB)
     +                        )
     +                     )
     +                  )
     +                -drho_dphi*phi(i,j)
     +                -drho_dpsi
     +                 *( dpsi_dphiA*phi(IA,JMax)
     +                   +dpsi_dphiB*phi(IB,JB)
     +                   +dpsi_dpsi0
     +                    *( dpsi0_dphi0*phi(I0,JMax)
     +                      +dpsi0_dphiA*phi(IA,JMax)
     +                      +dpsi0_dphiB*phi(IB,JB)
     +                      +dpsi0_dphim*phi(im,JMax)
     +                      +dpsi0_dphip*phi(ip,JMax)
     +                     )
     +                  )
     +               )
              b(k)= b(k)
     +             +4d0*pi
     +              *(-drho_dphi)
              bK0(k)=
     +             +4d0*pi
     +              *(-drho_dphiu
     +                 *( dphiu_dphi0
     +                   +dphiu_dC
     +                    *( dC_dphi0
     +                      +dC_dpsiA
     +                       *(+dpsiA_dpsi0A
     +                          *( dpsi0A_dphi0
     +                           )
     +                        )
     +                     )
     +                  )
     +                -drho_dC
     +                 *( dC_dphi0
     +                   +dC_dpsiA
     +                    *(+dpsiA_dpsi0A
     +                       *( dpsi0A_dphi0
     +                        )
     +                     )
     +                  )
     +                -drho_dpsi
     +                 *(+dpsi_dpsi0
     +                    *( dpsi0_dphi0
     +                     )
     +                  )
     +               )
              bKA(k)=
     +             +4d0*pi
     +              *(-drho_dphiu
     +                 *(+dphiu_dC
     +                    *(+dC_dphiA
     +                      +dC_dpsiA
     +                       *( dpsiA_dphiA
     +                         +dpsiA_dpsi0A
     +                          *(+dpsi0A_dphiA
     +                           )
     +                        )
     +                     )
     +                  )
     +                -drho_dC
     +                 *(+dC_dphiA
     +                   +dC_dpsiA
     +                    *( dpsiA_dphiA
     +                      +dpsiA_dpsi0A
     +                       *(+dpsi0A_dphiA
     +                        )
     +                     )
     +                  )
     +                -drho_dpsi
     +                 *( dpsi_dphiA
     +                   +dpsi_dpsi0
     +                    *(+dpsi0_dphiA
     +                     )
     +                  )
     +               )
              bKB(k)=
     +             +4d0*pi
     +              *(-drho_dphiu
     +                 *(+dphiu_dC
     +                    *(+dC_dpsiA
     +                       *(+dpsiA_dphiB
     +                         +dpsiA_dpsi0A
     +                          *(+dpsi0A_dphiB
     +                           )
     +                        )
     +                     )
     +                  )
     +                -drho_dC
     +                 *(+dC_dpsiA
     +                    *(+dpsiA_dphiB
     +                      +dpsiA_dpsi0A
     +                       *(+dpsi0A_dphiB
     +                        )
     +                     )
     +                  )
     +                -drho_dpsi
     +                 *(+dpsi_dphiB
     +                   +dpsi_dpsi0
     +                    *(+dpsi0_dphiB
     +                     )
     +                  )
     +               )
             bphim(k)=
     +             +4d0*pi
     +              *(-drho_dpsi
     +                 *(+dpsi_dpsi0
     +                    *(+dpsi0_dphim
     +                     )
     +                  )
     +               )
             kbphim(k)=kk(im,JMax)
             bphip(k)=
     +             +4d0*pi
     +              *(-drho_dpsi
     +                 *(+dpsi_dpsi0
     +                    *(+dpsi0_dphip
     +                     )
     +                  )
     +               )
	      kbphip(k)=kk(ip,JMax)
	    endif
          enddo
      enddo ! <== internal space ===
c     === forming of initial data for sparse matrix solver ==>
      NZ=0
      do k=1,KMax
          i=(k-1)/(JMax+1) 
	  j=k-1-i*(JMax+1)
          NZ=NZ+1
	  A(NZ)=b(k)
	  ir0(NZ)=k
	  ic0(NZ)=k
          if(i.ne.0) then 
              NZ=NZ+1
	      A(NZ)=bim(k)
	      ir0(NZ)=k
	      ic0(NZ)=kk(i-1,j)
          endif
          if(i.ne.IMax) then 
              NZ=NZ+1
	      A(NZ)=bip(k)
	      ir0(NZ)=k
	      ic0(NZ)=kk(i+1,j)
          endif
          if(j.ne.0) then 
              NZ=NZ+1
	      A(NZ)=bjm(k)
	      ir0(NZ)=k
	      ic0(NZ)=kk(i,j-1)
          endif
          if(j.ne.JMax) then 
              NZ=NZ+1
	      A(NZ)=bjp(k)
	      ir0(NZ)=k
	      ic0(NZ)=kk(i,j+1)
          endif
	  if((i.le.IA).and.(k.ge.kk(0,JMax))) then
              NZ=NZ+1
	      A(NZ)=bK0(k)
	      ir0(NZ)=k
	      ic0(NZ)=K0
              NZ=NZ+1
	      A(NZ)=bKA(k)
	      ir0(NZ)=k
	      ic0(NZ)=KA
              NZ=NZ+1
	      A(NZ)=bKB(k)
	      ir0(NZ)=k
	      ic0(NZ)=KB
              if(LogRapid) then
                  NZ=NZ+1
                  A(NZ)=bphim(k)
                  ir0(NZ)=k
                  ic0(NZ)=kbphim(k)
                  NZ=NZ+1
                  A(NZ)=bphip(k)
                  ir0(NZ)=k
                  ic0(NZ)=kbphip(k)
              endif
	  endif
      enddo
      do j=0,JMax
          NZ=NZ+1
          A(NZ)=b0(j)
	  ir0(NZ)=kk(0,JMax)
	  ic0(NZ)=kk(1,j)
      enddo
      do n=1,NZ
          if(LogInit) then
	      IRN(n)=ir0(n)
	      ICN(n)=ic0(n)
          else
	      IVECT(n)=ir0(n)
              JVECT(n)=ic0(n)
          endif
      enddo
      U=0.1d0
      LBLOCK=.TRUE.
      GROW=.TRUE.
      ABORT(1)=.TRUE.
      ABORT(2)=.TRUE.
      ABORT(3)=.FALSE.
      ABORT(4)=.FALSE.
      ! IFAIL=110
      IFAIL=010
      CALL X04ABF(1,NOUT)
      if(LogInit) then
***          CALL F01BRF(KMax,NZ,A,LICN,IRN,LIRN,ICN,U,IKEEP,IW,W,LBLOCK,
***     +                GROW,ABORT,IDISP,IFAIL)
          THR=1d-10
          CALL M28Y12(KMax,NZ,A,LICN,IRN,LIRN,ICN,U,IKEEP,IW,W,LBLOCK,
     +                GROW,ABORT,IDISP,IFAIL,THR)
          LogInit=.false.
c          LogInit=.true. ! only M28BYS (25 October)
c          NM28BYS=0
      else
          ETA=1d-4
          ! IFAIL=110
          IFAIL=010
***          CALL F01BSF(KMax,NZ,A,LICN,IVECT,JVECT,ICN,IKEEP,IW,W,GROW,
***     +                ETA,RPMIN,ABORT(4),IDISP,IFAIL)
          CALL M28BYS(KMax,NZ,A,LICN,IVECT,JVECT,ICN,IKEEP,IW,W,GROW,
     +                ETA,RPMIN,ABORT(4),IDISP,IFAIL)
c          NM28BYS=NM28BYS+1
c          if(NM28BYS.ge.2) LogInit=.true.
      endif
      MTYPE=1
      CALL M28CYS(KMax,A,LICN,ICN,IKEEP,ff,W,MTYPE,IDISP,RESID)
c     <== forming of initial data for sparse matrix solver ==>
      do k=1,KMax
          i=(k-1)/(JMax+1) 
	  j=k-1-i*(JMax+1)
          phi(i,j)=ff(k)
      enddo
      if(JB.eq.JMax) then ! thoroidal configurations
          ! let new  I0
          I00=I0
          potential_min=phi(I0,JMax)+Fpsi(degree(r(I0),delta))
          do i=IB+1,IA-1
              potential=phi(i,JMax)+Fpsi(degree(r(i),delta))
              if(potential.lt.potential_min) then
                  I0=i
                  potential_min=potential
              endif
          enddo
          if(I00.ne.I0) LogInit=.false.
      endif
      C=( phi(IA,JMax)+Fpsi(one)
     +   -phi(I0,JMax)*EOS(zero,3)/EOS(rho0,3)
     +  )
     +  /(one-EOS(zero,3)/EOS(rho0,3))
      phiu=EOS(rho0,3)/(C-phi(I0,JMax))
      return
      end

      integer Function kk(i,j)
      !*****************************************************************
      include 'statd.for'
      kk=1+i*(JMax+1)+j
      return
      end

      double precision Function Fpsi(om_)
      !*****************************************************************
      include 'accur.for'
      include 'statd.for'
      Common /phi/phi(0:IMax,0:JMax)
      Common /IB/IB/JB/JB

      Fpsi=(phi(IB,JB)-phi(IA,JMax))*Fpsi0(om_)
      return
      end

      double precision Function Fpsi0(om_)
      !*****************************************************************
      ! Fpsi0(r(I0))=zero
      ! Fpsi0(r(IA))-Fpsi0(degree(r(IB),delta)*sin(theta(JB)))=one
      include 'accur.for'
      include 'statd.for'
      Common /phi/phi(0:IMax,0:JMax)
      Common /delta/delta
      Common /IB/IB/JB/JB/I0/I0

      if(JB.eq.0) then
          iomb=0
      else
          iomb=IB
      endif
      if(om_.lt.one) then
          dr=one/IA
          im=degree(om_,one/delta)/dr
          ip=im+1
          omm=im*dr
          omp=ip*dr
          Fpsi0=( (omp-degree(om_,one/delta))/dr*phi(im,JMax)
     +           +(degree(om_,one/delta)-omm)/dr*phi(ip,JMax)
     +	         -phi(I0,JMax)
     +          )
     +          /(phi(IA,JMax)-phi(iomb,JMax))
      else
          Fpsi0=(phi(IA,JMax)-phi(I0,JMax))
     +          /(phi(IA,JMax)-phi(iomb,JMax))
      endif
      return
      end

      double precision Function Fdpsi_domega(om_)
      !*****************************************************************
      include 'accur.for'
      include 'statd.for'

      if(om_.lt.one) then
          dom=1d-5
          omm=max(zero,om_-half*dom)
          omp=min(one,om_+half*dom)
          Fdpsi_domega=(Fpsi(omp)-Fpsi(omm))/(omp-omm)
      else 
          Fdpsi_domega=zero
      endif
      return
      end

      double precision Function degree(x,alpha)
      !*****************************************************************
      include 'accur.for'

      if(alpha.eq.zero) then
          degree=one
      else 
          if(x.eq.zero) then
              degree=zero
          else
              degree=x**alpha
          endif
      endif
      return
      end

      Subroutine Integrals(CM,V,CJ,T,W,Eth,Pip,phi,r,theta,
     +                     C,phiu,rho0,omegaA)
      !*****************************************************************
      include 'accur.for'
      include 'statd.for'
      Dimension r(0:IMax),theta(0:JMax),
     +          phi(0:IMax,0:JMax)
      Common /pi/pi/IB/IB/JB/JB
      Common /delta/delta

      CM=zero
      V=zero
      CJ=zero
      T=zero
      W=zero
      Pip=zero
      Eth=zero
      do i=1,IA
          do j=1,JMax
              ra=half*(r(i-1)+r(i))
              thetaa=half*(theta(j-1)+theta(j))
              phia=(phi(i-1,j-1)+phi(i-1,j)+phi(i,j-1)+phi(i,j))/four
              psi=Fpsi(degree(ra,delta)*sin(thetaa))
              rho=EOS(phiu*(C-phia-psi),1)/rho0
     	      c1=-Fdpsi_domega(degree(ra,delta)*sin(thetaa))
     +            *degree(ra,delta)*sin(thetaa)
	      vphi=sqrt(max(zero,c1))
              dV=(degree(r(i),delta)**3-degree(r(i-1),delta)**3)/3d0
     +           *(cos(theta(j-1))-cos(theta(j)))
     +           *4d0*pi
	      CM=CM+rho*dV
	      if(rho.gt.zero) then
		  V=V+dV
	      endif
	      CJ= CJ
     +		 +rho*vphi
     +            *degree(ra,delta)*sin(thetaa)
     +            *dV
              p=EOS(phiu*(C-phia-psi),4)/(rho0*phiu)
              Pip=Pip+p*dV
	      T=T+rho*vphi**2/2d0*dV
	      W=W+rho*phia/2d0*dV
	      Eth=Eth+EOS(rho0*rho,5)/(rho0*phiu)*dV
              if((i.eq.IA).and.(j.eq.JMax)) then
                  omegaA=sqrt(max(zero,
     +                            c1/(degree(ra,delta)*sin(thetaa))
     +                           )
     +                       )
              endif
          enddo
      enddo
      return
      end
